% Computes the steady state of test.mod
% File sends "SS values" to twosector5v2.mod and updates the parameters
% that represent some steady state ratios (= "steady state parameters")

% Christoph Gortz

function [ys,check] = twosector5v7JT_steadystate(ys,exe)
global M_



%% 1. I load the values of the deep parameters in a loop.
%%
NumberOfParameters = M_.param_nbr;                            % Number of deep parameters.
for i = 1:NumberOfParameters                                  % Loop...
  paramname = deblank(M_.param_names(i,:));                   %    Get the name of parameter i. 
  eval([ paramname ' = M_.params(' int2str(i) ');']);         %    Get the value of parameter i.
end              


check = 0; % check=0 limits steady state values to real numbers (see  Dynare user guide page 71).



%% 2. Compute the values for the steady state parameters
%%
Lss = 1;
cbeta = 1/(1+constebeta/100);
ga = cga/100 ;
gv = cgv/100 ;
clandap = 1 + clandapaux;
clandaw = 1 + clandawaux;
%crho = -(cminrho + 1);
crho = 1/(cminrho - 1); % restricting cminrho in (0,1)
%cga = constepinfi - constepinfc -  (alfac-1)/(1-alfai)*cgv;
%ga = cga/100;


% For start values in vector x0:
% Steady state values without adjustment costs a la Huffman and Wynne
% (used for initial guess)

iX1aux = (1-alfai)^(alfai-1)*alfai^-alfai*(1-alfac)^(1-alfai)*alfac^(alfac*(1-alfai)/(1-alfac));
% steady state rental rate
icrk = (( exp((1/(1-alfai))*gv)/cbeta - (1-cdeltai) )*iX1aux*(1/clandap)^((alfac-alfai)/(1-alfac)) )^((1-alfac)/(1-alfai)); 
% steady state relative price of investment
icpi = iX1aux* icrk^((alfai-alfac)/(1-alfac))*(1/clandap)^((alfac-alfai)/(1-alfac));
% steady state wage
icw = ( (1/clandap)*(1-alfac)^(1-alfac)*alfac^alfac*icrk^-alfac )^(1/(1-alfac));
% steady state capital labor ratio in consumption sector
iKcLc = icw/icrk*(alfac/(1-alfac)); 
% steady state capital labor ratio in investment sector
iKiLi = icw/icrk*(alfai/(1-alfai)); 
% Labour
iLc = Lss*( 1+ nvalueic*(iKcLc^alfac)/(iKiLi^alfai)*icpi^-1 )^-1;
iLi = Lss - iLc;
% steady state capital in consumption sector
iKc = iKcLc*iLc;
% steady state capital in investment sector
iKi = iKiLi*iLi;
% steady state investment sector output
iiss = (1/clandap)*iKiLi^alfai*iLi; 
% steady state investement in I sector:
iiiss = iKi*(exp(1/(1-alfai)*gv) - (1-cdeltai));
% steady state consumption in I sector:
iicss = iKc*(exp(1/(1-alfai)*gv) - (1-cdeltac));

% Calculate values of "steady state parameters"
% Steady state equations for SStwosectorv5.mod (2 Sector model with
% adjustment costs a la Huffman and Wynne)

% Relabeling the variables:
% ki = x(1), kc = x(2), ii = x(3), ic = x(4) rik = x(5) rck = x(6) Li = x(7) Lc = x(8) w = x(9)  pi = x(10)

% % F contains the model equations in SS (all equal to zero):
 F = @(x) [
% % Equation ki
 x(1) - x(9)/x(5) * alfai/(1-alfai) * x(7);
% % Equation kc
 x(2) - x(9)/x(6) * alfac/(1-alfac) * x(8);
% % Equation w 
 x(9) - ( 1/clandap * (1-alfac)^(1-alfac)*alfac^alfac*x(6)^(-alfac) )^(1/(1-alfac));
% % Equation pi
%x(10) - ( 1/clandapaux *(1-alfac)^(1-alfac)*alfac^alfac*x(6)^(-alfac) )^((1-alfai)/(1-alfac)) / (1/clandapaux * (1-alfai)^(1-alfai)*alfai^alfai*(x(5))^(-alfai) );
%x(10) - ( (1-alfai)^(alfai-1)*alfai^(-alfai) )/( (1-alfac)^(alfac-1)*alfac^(-alfac) ) *(x(9))^(alfac-alfai) * (x(6))^(-alfac)*(x(5))^alfai;
x(10) - ( 1/clandap*(1-alfac)^(1-alfac)*alfac^alfac * ( (exp(1/(1-alfai)*gv))/cbeta - (1-cdeltac) )^-alfac * ((x(3)^(-crho)+x(4)^(-crho))^(-1/crho -1)*x(4)^(-crho-1))^-alfac )  /  (  1/clandap*(1-alfai)^(1-alfai)*alfai^alfai * ( (exp(1/(1-alfai)*gv))/cbeta - (1-cdeltai) )^-alfai * ((x(3)^(-crho)+x(4)^(-crho))^(-1/crho -1)*x(3)^(-crho-1))^-alfai    )^((1-alfac)/(1-alfai));
% % Equation Lc
 x(8) - (1+nvalueic*(x(2)/x(8))^alfac/((x(1)/x(7))^alfai) * x(10)^(-1) )^(-1);
% % Equation Li
 x(7) - Lss + x(8);
% % Equation rki
 x(5) - ( (exp(1/(1-alfai)*gv))/cbeta - (1-cdeltai) )*x(10) * (x(3)^(-crho)+x(4)^(-crho))^(-1/crho -1)*x(3)^(-crho-1);
% % Equation rkc
 x(6) - ( (exp(1/(1-alfai)*gv))/cbeta - (1-cdeltac) )*x(10) * (x(3)^(-crho)+x(4)^(-crho))^(-1/crho -1)*x(4)^(-crho-1);
% % Equation ii
 x(3) - x(1)*(1-exp(-1/(1-alfai)*gv)*(1-cdeltai));
% % Equation ic
 x(4) - x(2)*(1-exp(-1/(1-alfai)*gv)*(1-cdeltac));
 ];

 x0 = [iKi; iKc; iiiss; iicss; icrk; icrk; iLi; iLc; icw; icpi;];     % Initial guess for start values
 options = optimset('Display','off','TolFun',1e-12,'TolX',1e-12);
 x = fsolve(F,x0,options);  % x contains the variables x(1)-x(10)

X1aux = (1-alfai)^(alfai-1)*alfai^-alfai*(1-alfac)^(1-alfai)*alfac^(alfac*(1-alfai)/(1-alfac));
crki = x(5);
crkc = x(6);
cpi = x(10);
cw = x(9);
%KcLc = x(8)*x(2);
KcLc = x(2)/x(8);
%KiLi = x(7)*x(1);
KiLi = x(1)/x(7);
Lc = x(8);
Li = x(7);
css = (1/clandap)*KcLc^alfac*Lc;
Kc = x(2);
Ki = x(1);
K = x(2) + x(1);
Kibar = Ki*exp((1/(1-alfai))*gv);
Kcbar = Kc*exp((1/(1-alfai))*gv);
iiss = x(3);
icss = x(4);
iss = (1/clandap)*KiLi^alfai*Li; 
yss = (css+cpi*iss);

% Steady state of price of capital in I sector
qiss = (iiss^(-crho) + icss^(-crho))^(-1/crho-1) * iiss^(-crho-1) * cpi;
% Steady state of price of capital in C sector
qcss = (iiss^(-crho) + icss^(-crho))^(-1/crho-1) * icss^(-crho-1) * cpi;

rbiss = (crki + qiss*(1-cdeltai))/qiss;
rbcss = (crkc + qcss*(1-cdeltai))/qcss;

%% 5. Assign SS variables for vector ys
%%
 
labis = 0;
labcs = 0;
labs =0;
kcs = 0;
kis = 0;
kpis = 0;
kpcs = 0;
mccs = 0;
mcis = 0;
zcapis = 0;
zcapcs = 0;
rkis = 0;
rkcs = 0;
pis = 0;
cs = 0;
lams = 0;
phifis = 0;
phifcs = 0;
inveis = 0;
invecs = 0;
inves = 0;
ysss = 0;
pinfcs = 0;
pinfis = 0;
ws =0;
rs = 0;
zs = 0; 
vs = 0; 
z_ls = 0;
v_ls = 0;
bs = 0;
gs = 0;
qss = 0;
mss = 0; 
sws = 0;
ewmas = 0;
dys = 0 + 0 + alfac/(1-alfai)*0 + (cga + alfac/(1-alfai)*cgv)*0;  % Observables...
dcs = 0-0 + 0 + alfac/(1-alfai)*0 + (cga + alfac/(1-alfai)*cgv)*0;
dinves = 0-0 + 1/(1-alfai)*0 + ( 1/(1-alfai)*cgv)*0;
dpis = 0 + 0 + (alfac-1)/(1-alfai)*0 + (cga+ (alfac-1)/(1-alfai)*cgv)*0;
dws = 0 + 0 + alfac/(1-alfai)*0 + (cga + alfac/(1-alfai)*cgv)*0;
labobss = 0;
robss = 0 +( log(constepinfc/100 + 1) - log(cbeta) + (cga + alfac/(1-alfai)*cgv))*0;
pinfobscs = 1*(0) + constepinfc*0;
pinfobsis = 1*(0) + constepinfi*0 ;
%pinfobsis = pinfobscs + dpis;
rrates = 0;     % ... end observables.
epinfmacs = 0;    
epinfmais = 0;
spinfcs = 0;
spinfis = 0;
zcapifs = zcapis;    % flexible economy variables. Should have same ss as above.
zcapcfs = zcapcs;
rkcfs = rkcs;
rkifs = rkis;
kpifs = kpis;
kpcfs = kpcs;
kifs = kis;
kcfs = kcs;
lamfs = lams;
phififs = phifis;
phifcfs = phifcs;
cfs = cs;
invefs = inves;
inveifs = inveis;
invecfs = invecs;
yfs = ysss;
labfs = labs;
labifs = labis;
labcfs = labcs;
pifs = pis;
wfs = ws;
spis = 0;
srps = 0;
spkcs = 0;
spkis = 0;
epkc1aux = 0;
epkc4aux = 0;
epkc8aux = 0;
epki1aux = 0;
epki4aux = 0;
epki8aux = 0;
ez1aux = 0;
ez4aux = 0;
ez8aux = 0;
ev1aux = 0;
ev4aux = 0;
ev8aux = 0;
rbcs = 0;
rbis = 0;
qcs = 0;
qis = 0;
etfp4aux = 0; 
etfp8aux = 0;
tfps = 0;


ys = [ labis; labcs; labs; kcs; kis; kpis; kpcs; mccs; mcis; zcapis; zcapcs; rkis; rkcs;
pis; lams; phifis; phifcs; cs; inveis; invecs; inves; ysss; pinfcs; pinfis; ws; rs; zs;
vs; z_ls; v_ls; bs; gs; qss; mss; sws; ewmas; dys; dcs; dinves;    dpis;     dws; labobss; robss; pinfobscs; 
pinfobsis; rrates; epinfmacs; epinfmais; spinfcs; spinfis;  zcapifs; zcapcfs; 
rkcfs; rkifs; kpifs; kpcfs; kifs; kcfs; lamfs; phififs; phifcfs; cfs; invefs; inveifs;
invecfs; yfs; labfs; labifs; labcfs; pifs; wfs; spis; srps;
spkcs; spkis; epkc1aux; epkc4aux; epkc8aux; epki1aux; epki4aux; epki8aux;
ez1aux; ez4aux; ez8aux; ev1aux; ev4aux; ev8aux; rbcs; rbis; qcs; qis; tfps; etfp4aux; etfp8aux;
];

%% 6. update the "steady state parameters" in the global workspace  
%% steady state parameters: X1aux crki crkc cpi cw KcLc KiLi Lc Li css iss
%% Kc Ki K Kibar Kcbar yss iiss icss

id1 = strmatch('X1aux',M_.param_names,'exact');     % Get the index of crki in M_.params
id2 = strmatch('crki',M_.param_names,'exact');      % Get the index of crkc in M_.params
id3 = strmatch('crkc',M_.param_names,'exact');                   
id4 = strmatch('cpi',M_.param_names,'exact');
id5 = strmatch('cw',M_.param_names,'exact');                   
id6 = strmatch('KcLc',M_.param_names,'exact');
id7 = strmatch('KiLi',M_.param_names,'exact');
id8 = strmatch('Lc',M_.param_names,'exact');
id9 = strmatch('Li',M_.param_names,'exact');
id10 = strmatch('css',M_.param_names,'exact');
id11 = strmatch('iss',M_.param_names,'exact');         
id12 = strmatch('Kc',M_.param_names,'exact');           
id13 = strmatch('Ki',M_.param_names,'exact');                   
id14 = strmatch('K',M_.param_names,'exact');
id15 = strmatch('Kibar',M_.param_names,'exact');                   
id16 = strmatch('Kcbar',M_.param_names,'exact');
id17 = strmatch('yss',M_.param_names,'exact');
id18 = strmatch('iiss',M_.param_names,'exact');
id19 = strmatch('icss',M_.param_names,'exact');
id20 = strmatch('qcss',M_.param_names,'exact');         
id21 = strmatch('qiss',M_.param_names,'exact');    
M_.params(id1) = X1aux;
M_.params(id2) = crki;
M_.params(id3) = crkc;
M_.params(id4) = cpi;
M_.params(id5) = cw;
M_.params(id6) = KcLc;
M_.params(id7) = KiLi;
M_.params(id8) = Lc;
M_.params(id9) = Li;
M_.params(id10) = css;
M_.params(id11) = iss;
M_.params(id12) = Kc;
M_.params(id13) = Ki;
M_.params(id14) = K;
M_.params(id15) = Kibar;
M_.params(id16) = Kcbar;
M_.params(id17) = yss;
M_.params(id18) = iiss;
M_.params(id19) = icss;
M_.params(id20) = qcss;
M_.params(id21) = qiss;
